knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_and_barriers', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan', # todo: do we need this??
    'barriers_self_cw', 
    'facilitators_self_cw', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'barriers_self_cb', 
    'facilitators_self_cb', 
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own actionplan",
    'Partner actionplan',
    "Daily barriers",
    "Daily facilitators",
    "Daily wear time",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    'Mean barriers',
    "Mean facilitators",
    "Mean wear time"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,13),
  "Between-Person Effects" = c(14,22),
  "Random Effects" = c(23, 29), 
  "Additional Parameters" = c(30,30)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,13+5),
  "Between-Person Effects" = c(14+5,22+5),
  "Random Effects" = c(23+5, 29+5), 
  "Additional Parameters" = c(30+5,30+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.33 [1.29, 1.38]         1.15      0.75
##  persuasion_partner_cw 1.08 [1.06, 1.13]         1.04      0.92
##       pressure_self_cw 1.05 [1.02, 1.09]         1.02      0.96
##    pressure_partner_cw 1.06 [1.04, 1.10]         1.03      0.94
##        pushing_self_cw 1.03 [1.01, 1.08]         1.02      0.97
##     pushing_partner_cw 1.14 [1.11, 1.19]         1.07      0.87
##     persuasion_self_cb 1.04 [1.02, 1.09]         1.02      0.96
##  persuasion_partner_cb 3.82 [3.63, 4.02]         1.95      0.26
##       pressure_self_cb 3.78 [3.59, 3.97]         1.94      0.26
##    pressure_partner_cb 2.40 [2.29, 2.51]         1.55      0.42
##        pushing_self_cb 2.36 [2.26, 2.47]         1.54      0.42
##     pushing_partner_cb 3.37 [3.21, 3.54]         1.83      0.30
##       barriers_self_cw 3.37 [3.21, 3.54]         1.83      0.30
##       barriers_self_cb 1.26 [1.22, 1.31]         1.12      0.79
##   facilitators_self_cw 1.12 [1.09, 1.16]         1.06      0.89
##   facilitators_self_cb 1.27 [1.23, 1.32]         1.13      0.79
##                    day 1.08 [1.06, 1.13]         1.04      0.92
##  Tolerance 95% CI
##      [0.72, 0.78]
##      [0.89, 0.95]
##      [0.92, 0.98]
##      [0.91, 0.96]
##      [0.92, 0.99]
##      [0.84, 0.90]
##      [0.92, 0.98]
##      [0.25, 0.28]
##      [0.25, 0.28]
##      [0.40, 0.44]
##      [0.40, 0.44]
##      [0.28, 0.31]
##      [0.28, 0.31]
##      [0.76, 0.82]
##      [0.86, 0.92]
##      [0.76, 0.81]
##      [0.89, 0.95]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         91%
##     lognormal          9%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 7 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 7, observations = 1776, p-value =
## 0.04
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003378378
## sample estimates:
## outlier frequency (expected: 0.00162725225225225 ) 
##                                        0.003941441
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report,
## rope_range = rope_range_continuous, : Coefficients were exponentiated.
## Double check if this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 2.46*** 0.61 [ 1.48, 4.02] 1.000 [0.84, 1.20] 0.003 22.432 Strong Evidence 1.000 34292 30616 49.83*** 3.96 [42.49, 58.35] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 12383 20018
Within-Person Effects
Daily individual’s experienced persuasion 1.46*** 0.14 [ 1.22, 1.78] 1.000 [0.84, 1.20] 0.015 >100 Overwhelming Evidence 1.000 54434 30868 1.02 0.03 [ 0.97, 1.09] 0.814 [0.93, 1.08] 0.961 0.016 Very Strong Evidence for Null 1.000 28321 32412
Daily partner’s experienced persuasion 1.32** 0.13 [ 1.10, 1.65] 0.999 [0.84, 1.20] 0.145 4.143 Moderate Evidence 1.000 44672 28708 1.02 0.02 [ 0.98, 1.07] 0.841 [0.93, 1.08] 0.986 0.015 Very Strong Evidence for Null 1.000 38484 32300
Daily individual’s experienced pressure 0.80 0.20 [ 0.47, 1.46] 0.810 [0.84, 1.20] 0.353 0.208 Moderate Evidence for Null 1.000 34482 23563 0.88 0.06 [ 0.77, 1.01] 0.969 [0.93, 1.08] 0.208 0.076 Strong Evidence for Null 1.000 36562 29427
Daily partner’s experienced pressure 1.50 0.58 [ 0.81, 4.82] 0.896 [0.84, 1.20] 0.225 0.313 Weak Evidence for Null 1.000 26034 22753 0.95 0.05 [ 0.86, 1.04] 0.851 [0.93, 1.08] 0.716 0.014 Very Strong Evidence for Null 1.000 44427 33916
Daily individual’s experienced pushing 1.15 0.17 [ 0.86, 1.57] 0.838 [0.84, 1.20] 0.586 0.118 Moderate Evidence for Null 1.000 40066 29999 0.97 0.03 [ 0.91, 1.03] 0.837 [0.93, 1.08] 0.924 0.013 Very Strong Evidence for Null 1.000 40647 37344
Daily partner’s experienced pushing 1.58** 0.27 [ 1.17, 2.43] 0.998 [0.84, 1.20] 0.035 6.351 Moderate Evidence 1.000 32077 24162 0.97 0.03 [ 0.91, 1.02] 0.875 [0.93, 1.08] 0.920 0.016 Very Strong Evidence for Null 1.000 48349 33082
Day 1.00 0.22 [ 0.65, 1.54] 0.501 [0.84, 1.20] 0.591 0.109 Moderate Evidence for Null 1.000 72747 30791 1.00 0.07 [ 0.87, 1.14] 0.516 [0.93, 1.08] 0.737 0.008 Very Strong Evidence for Null 1.000 78601 30176
Own actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.96 0.09 [ 0.80, 1.15] 0.680 [0.84, 1.20] 0.919 0.053 Strong Evidence for Null 1.000 66692 32610 1.21*** 0.04 [ 1.13, 1.30] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 80200 30400
Daily facilitators 1.40*** 0.13 [ 1.17, 1.68] 1.000 [0.84, 1.20] 0.046 28.747 Strong Evidence 1.000 70661 33120 1.20*** 0.03 [ 1.14, 1.26] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 86769 28603
Daily wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.86 0.47 [ 0.28, 2.49] 0.611 [0.84, 1.20] 0.248 0.282 Moderate Evidence for Null 1.000 23398 28030 1.07 0.17 [ 0.78, 1.46] 0.655 [0.93, 1.08] 0.337 0.016 Very Strong Evidence for Null 1.000 23037 30106
Mean partner’s experienced persuasion 1.31 0.72 [ 0.43, 3.81] 0.690 [0.84, 1.20] 0.223 0.313 Weak Evidence for Null 1.000 23206 27678 0.92 0.15 [ 0.67, 1.26] 0.696 [0.93, 1.08] 0.325 0.018 Very Strong Evidence for Null 1.000 22566 28354
Mean individual’s experienced pressure 0.31 0.23 [ 0.07, 1.27] 0.947 [0.84, 1.20] 0.056 1.308 Weak Evidence 1.000 34185 32075 0.64 0.25 [ 0.29, 1.40] 0.869 [0.93, 1.08] 0.081 0.045 Strong Evidence for Null 1.000 9237 16406
Mean partner’s experienced pressure 0.25 0.19 [ 0.05, 1.05] 0.971 [0.84, 1.20] 0.035 2.069 Weak Evidence 1.000 32700 30664 0.49 0.19 [ 0.23, 1.08] 0.962 [0.93, 1.08] 0.031 0.112 Moderate Evidence for Null 1.000 9095 15701
Mean individual’s experienced pushing 1.13 0.91 [ 0.23, 5.71] 0.560 [0.84, 1.20] 0.178 0.409 Weak Evidence for Null 1.000 29830 31903 1.49 0.43 [ 0.83, 2.64] 0.911 [0.93, 1.08] 0.082 0.043 Strong Evidence for Null 1.000 12430 20999
Mean partner’s experienced pushing 1.42 1.16 [ 0.28, 7.17] 0.667 [0.84, 1.20] 0.159 0.453 Weak Evidence for Null 1.000 29928 31218 1.70 0.50 [ 0.94, 3.03] 0.961 [0.93, 1.08] 0.043 0.084 Strong Evidence for Null 1.000 12609 21545
Mean barriers 0.57* 0.13 [ 0.37, 0.90] 0.992 [0.84, 1.20] 0.049 2.229 Weak Evidence 1.000 78953 29599 0.99 0.07 [ 0.86, 1.15] 0.535 [0.93, 1.08] 0.693 0.010 Very Strong Evidence for Null 1.000 53450 32820
Mean facilitators 1.21 0.19 [ 0.89, 1.64] 0.886 [0.84, 1.20] 0.468 0.161 Moderate Evidence for Null 1.000 71293 31096 1.11* 0.05 [ 1.01, 1.22] 0.985 [0.93, 1.08] 0.290 0.127 Moderate Evidence for Null 1.000 40796 32579
Mean wear time NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.25 0.19 [0.94, 1.69] NA NA NA NA NA 1.000 26033 28431 0.30 0.04 [0.22, 0.40] NA NA NA NA NA 1.000 16013 25674
sd(Daily individual’s experienced persuasion) 0.16 0.12 [0.01, 0.44] NA NA NA NA NA 1.000 18981 21658 0.11 0.03 [0.06, 0.17] NA NA NA NA NA 1.000 20545 19744
sd(Daily partner’s experienced persuasion) 0.24 0.14 [0.02, 0.53] NA NA NA NA NA 1.000 14487 16455 0.06 0.03 [0.01, 0.12] NA NA NA NA NA 1.000 15927 14482
sd(Daily individual’s experienced pressure) 0.43 0.39 [0.02, 1.79] NA NA NA NA NA 1.000 14462 21703 0.09 0.08 [0.00, 0.34] NA NA NA NA NA 1.000 16675 22880
sd(Daily partner’s experienced pressure) 0.78 0.60 [0.04, 2.35] NA NA NA NA NA 1.000 14322 19555 0.05 0.05 [0.00, 0.19] NA NA NA NA NA 1.000 24187 23835
sd(Daily individual’s experienced pushing) 0.45 0.23 [0.05, 0.96] NA NA NA NA NA 1.000 12397 13902 0.06 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 15318 18397
sd(Daily partner’s experienced pushing) 0.42 0.29 [0.03, 1.14] NA NA NA NA NA 1.000 11394 18073 0.05 0.04 [0.00, 0.13] NA NA NA NA NA 1.000 16710 18837
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.66 0.01 [0.63, 0.69] NA NA NA NA NA 1.000 66013 28624 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.76), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.8). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...
## Extracting .afm files from .ttf files...
## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 131 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 204 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0077

$persuasion_partner_cw

## Warning: Removed 101 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 198 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.00703

$pressure_self_cw

## Warning: Removed 65 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Warning: Removed 91 rows containing missing values or values outside the scale
## range (`geom_line()`).
## Picking joint bandwidth of 0.0126

$pressure_partner_cw

## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.028

$pushing_self_cw

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00906

$pushing_partner_cw

## Picking joint bandwidth of 0.013

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0     -0.1      0.08    -0.23     0.03       8.25
##   Post.Prob Star
## 1      0.89     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner + 
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.12 [1.08, 1.18]         1.06      0.89
##  persuasion_partner_cw 1.16 [1.11, 1.22]         1.08      0.86
##       pressure_self_cw 1.08 [1.04, 1.14]         1.04      0.93
##    pressure_partner_cw 1.11 [1.07, 1.17]         1.05      0.90
##        pushing_self_cw 1.17 [1.13, 1.23]         1.08      0.85
##     pushing_partner_cw 1.16 [1.12, 1.22]         1.08      0.86
##       pressure_self_cb 3.55 [3.32, 3.80]         1.88      0.28
##    pressure_partner_cb 3.70 [3.46, 3.96]         1.92      0.27
##       barriers_self_cw 1.34 [1.28, 1.42]         1.16      0.75
##       barriers_self_cb 1.16 [1.12, 1.22]         1.08      0.86
##   facilitators_self_cw 1.34 [1.28, 1.41]         1.16      0.75
##   facilitators_self_cb 1.13 [1.09, 1.19]         1.06      0.88
##                    day 1.07 [1.03, 1.13]         1.03      0.94
##       weartime_self_cw 1.06 [1.03, 1.12]         1.03      0.94
##       weartime_self_cb 1.22 [1.17, 1.28]         1.10      0.82
##  Tolerance 95% CI
##      [0.85, 0.93]
##      [0.82, 0.90]
##      [0.88, 0.96]
##      [0.86, 0.94]
##      [0.81, 0.89]
##      [0.82, 0.89]
##      [0.26, 0.30]
##      [0.25, 0.29]
##      [0.71, 0.78]
##      [0.82, 0.89]
##      [0.71, 0.78]
##      [0.84, 0.92]
##      [0.89, 0.97]
##      [0.89, 0.97]
##      [0.78, 0.85]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 6.70 [6.23, 7.21]         2.59      0.15
##  persuasion_partner_cb 7.21 [6.71, 7.76]         2.69      0.14
##        pushing_self_cb 5.74 [5.35, 6.17]         2.40      0.17
##     pushing_partner_cb 5.68 [5.29, 6.10]         2.38      0.18
##  Tolerance 95% CI
##      [0.14, 0.16]
##      [0.13, 0.15]
##      [0.16, 0.19]
##      [0.16, 0.19]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 13, observations = 1594, p-value =
## 0.02
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0006273526 0.0062735257
## sample estimates:
## outlier frequency (expected: 0.00308657465495609 ) 
##                                        0.008155583
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 122.48*** 7.80 [108.02, 139.14] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 9877 18232
Within-Person Effects
Daily individual’s experienced persuasion 1.02 0.02 [ 0.98, 1.06] 0.864 [0.94, 1.06] 0.991 0.013 Very Strong Evidence for Null 1.000 36109 33533
Daily partner’s experienced persuasion 1.03 0.02 [ 0.99, 1.07] 0.945 [0.94, 1.06] 0.969 0.025 Very Strong Evidence for Null 1.000 44766 32167
Daily individual’s experienced pressure 0.95 0.05 [ 0.86, 1.05] 0.853 [0.94, 1.06] 0.582 0.013 Very Strong Evidence for Null 1.000 43033 30808
Daily partner’s experienced pressure 0.96 0.04 [ 0.90, 1.04] 0.833 [0.94, 1.06] 0.747 0.010 Very Strong Evidence for Null 1.000 57513 31600
Daily individual’s experienced pushing 1.02 0.03 [ 0.97, 1.08] 0.807 [0.94, 1.06] 0.905 0.011 Very Strong Evidence for Null 1.000 41918 33775
Daily partner’s experienced pushing 1.01 0.02 [ 0.96, 1.05] 0.618 [0.94, 1.06] 0.987 0.007 Very Strong Evidence for Null 1.000 56068 33129
Day 0.97 0.05 [ 0.88, 1.07] 0.744 [0.94, 1.06] 0.701 0.007 Very Strong Evidence for Null 1.000 78736 29523
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.01 0.02 [ 0.96, 1.05] 0.589 [0.94, 1.06] 0.994 0.006 Very Strong Evidence for Null 1.000 81630 29963
Daily facilitators 1.06** 0.02 [ 1.02, 1.10] 0.999 [0.94, 1.06] 0.576 0.979 Weak Evidence for Null 1.000 81119 29545
Daily wear time 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.06] 1.000 1.635 Weak Evidence 1.000 83564 29235
Between-Person Effects
Mean individual’s experienced persuasion 1.02 0.17 [ 0.73, 1.42] 0.551 [0.94, 1.06] 0.293 0.015 Very Strong Evidence for Null 1.000 9540 16159
Mean partner’s experienced persuasion 0.89 0.15 [ 0.64, 1.25] 0.752 [0.94, 1.06] 0.230 0.020 Very Strong Evidence for Null 1.000 9520 15991
Mean individual’s experienced pressure 1.07 0.20 [ 0.74, 1.56] 0.645 [0.94, 1.06] 0.244 0.011 Very Strong Evidence for Null 1.000 12821 23088
Mean partner’s experienced pressure 1.02 0.19 [ 0.70, 1.45] 0.536 [0.94, 1.06] 0.266 0.011 Very Strong Evidence for Null 1.000 11811 19751
Mean individual’s experienced pushing 1.11 0.27 [ 0.68, 1.80] 0.666 [0.94, 1.06] 0.182 0.014 Very Strong Evidence for Null 1.000 12262 22240
Mean partner’s experienced pushing 1.29 0.31 [ 0.80, 2.08] 0.855 [0.94, 1.06] 0.116 0.024 Very Strong Evidence for Null 1.001 11897 20336
Mean barriers 0.81*** 0.04 [ 0.72, 0.89] 1.000 [0.94, 1.06] 0.002 10.406 Strong Evidence 1.000 64174 29854
Mean facilitators 0.95 0.03 [ 0.89, 1.02] 0.906 [0.94, 1.06] 0.650 0.023 Very Strong Evidence for Null 1.000 44618 32842
Mean wear time 1.00* 0.00 [ 1.00, 1.00] 0.980 [0.94, 1.06] 1.000 0.080 Strong Evidence for Null 1.000 39933 32096
Random Effects
sd(Intercept) 0.34 0.05 [0.26, 0.45] NA NA NA NA NA 1.000 12542 21650
sd(Daily individual’s experienced persuasion) 0.05 0.02 [0.02, 0.09] NA NA NA NA NA 1.000 17237 12946
sd(Daily partner’s experienced persuasion) 0.05 0.02 [0.00, 0.10] NA NA NA NA NA 1.000 11835 13979
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.23] NA NA NA NA NA 1.000 15112 21630
sd(Daily partner’s experienced pressure) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 25444 22855
sd(Daily individual’s experienced pushing) 0.09 0.03 [0.02, 0.16] NA NA NA NA NA 1.000 11352 8267
sd(Daily partner’s experienced pushing) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 15057 19346
Additional Parameters
sigma 0.54 0.01 [0.52, 0.56] NA NA NA NA NA 1.000 61564 28538 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.84), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.78). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.19 [1.15,  1.25]         1.09      0.84
##  persuasion_partner_cw 1.13 [1.09,  1.19]         1.07      0.88
##       pressure_self_cw 1.13 [1.09,  1.19]         1.07      0.88
##    pressure_partner_cw 1.07 [1.04,  1.13]         1.03      0.93
##        pushing_self_cw 1.27 [1.22,  1.33]         1.13      0.79
##     pushing_partner_cw 1.13 [1.09,  1.19]         1.06      0.88
##       pressure_self_cb 4.82 [4.51,  5.17]         2.20      0.21
##    pressure_partner_cb 4.76 [4.45,  5.10]         2.18      0.21
##       barriers_self_cw 1.34 [1.28,  1.41]         1.16      0.75
##       barriers_self_cb 1.13 [1.09,  1.18]         1.06      0.89
##   facilitators_self_cw 1.33 [1.27,  1.40]         1.15      0.75
##   facilitators_self_cb 1.09 [1.05,  1.14]         1.04      0.92
##                    day 1.07 [1.04,  1.13]         1.04      0.93
##  Tolerance 95% CI
##      [0.80, 0.87]
##      [0.84, 0.91]
##      [0.84, 0.91]
##      [0.88, 0.96]
##      [0.75, 0.82]
##      [0.84, 0.92]
##      [0.19, 0.22]
##      [0.20, 0.22]
##      [0.71, 0.78]
##      [0.84, 0.92]
##      [0.72, 0.79]
##      [0.87, 0.95]
##      [0.88, 0.96]
## 
## Moderate Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 9.80 [9.11, 10.54]         3.13      0.10
##  persuasion_partner_cb 9.86 [9.17, 10.60]         3.14      0.10
##        pushing_self_cb 7.30 [6.80,  7.84]         2.70      0.14
##     pushing_partner_cb 7.35 [6.85,  7.90]         2.71      0.14
##  Tolerance 95% CI
##      [0.09, 0.11]
##      [0.09, 0.11]
##      [0.13, 0.15]
##      [0.13, 0.15]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         44%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 11, observations = 1776, p-value =
## 0.00112
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.003095804 0.011055146
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.006193694
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.80*** 0.11 [ 3.57, 4.01] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 6797 12185
Within-Person Effects
Daily individual’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.631 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 49498 28462
Daily partner’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.687 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 55442 28931
Daily individual’s experienced pressure -0.04 0.06 [-0.18, 0.08] 0.751 [-0.11, 0.11] 0.841 0.006 Very Strong Evidence for Null 1.000 44958 27175
Daily partner’s experienced pressure -0.06 0.06 [-0.19, 0.06] 0.845 [-0.11, 0.11] 0.789 0.009 Very Strong Evidence for Null 1.000 41645 27925
Daily individual’s experienced pushing 0.04 0.04 [-0.04, 0.11] 0.841 [-0.11, 0.11] 0.973 0.008 Very Strong Evidence for Null 1.000 50741 33820
Daily partner’s experienced pushing 0.11* 0.04 [ 0.03, 0.18] 0.993 [-0.11, 0.11] 0.550 0.156 Moderate Evidence for Null 1.000 36214 28776
Day 0.24** 0.08 [ 0.09, 0.39] 0.999 [-0.11, 0.11] 0.048 0.602 Weak Evidence for Null 1.000 66580 29222
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers -0.14*** 0.03 [-0.20, -0.07] 1.000 [-0.11, 0.11] 0.203 7.138 Moderate Evidence 1.000 79478 29981
Daily facilitators 0.12*** 0.03 [ 0.06, 0.17] 1.000 [-0.11, 0.11] 0.382 19.986 Strong Evidence 1.000 73684 29956
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.14 0.29 [-0.44, 0.73] 0.689 [-0.11, 0.11] 0.263 0.015 Very Strong Evidence for Null 1.000 6587 13866
Mean partner’s experienced persuasion 0.45 0.29 [-0.13, 1.03] 0.935 [-0.11, 0.11] 0.095 0.047 Strong Evidence for Null 1.000 6385 13223
Mean individual’s experienced pressure -0.16 0.30 [-0.76, 0.43] 0.707 [-0.11, 0.11] 0.251 0.010 Very Strong Evidence for Null 1.000 8589 17922
Mean partner’s experienced pressure -0.32 0.30 [-0.91, 0.27] 0.858 [-0.11, 0.11] 0.166 0.016 Very Strong Evidence for Null 1.000 8134 16845
Mean individual’s experienced pushing 0.27 0.42 [-0.58, 1.11] 0.740 [-0.11, 0.11] 0.167 0.015 Very Strong Evidence for Null 1.001 7935 14546
Mean partner’s experienced pushing -0.09 0.42 [-0.93, 0.76] 0.586 [-0.11, 0.11] 0.203 0.012 Very Strong Evidence for Null 1.001 7770 14300
Mean barriers -0.48*** 0.08 [-0.64, -0.32] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 63689 32054
Mean facilitators 0.29*** 0.06 [ 0.18, 0.40] 1.000 [-0.11, 0.11] 0.001 >100 Overwhelming Evidence 1.000 36704 30529
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.62 0.08 [0.49, 0.81] NA NA NA NA NA 1.000 9809 18127
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 13113 16901
sd(Daily partner’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 14690 18195
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.26] NA NA NA NA NA 1.000 20756 20625
sd(Daily partner’s experienced pressure) 0.09 0.07 [0.00, 0.29] NA NA NA NA NA 1.000 17910 17916
sd(Daily individual’s experienced pushing) 0.08 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 12798 15220
sd(Daily partner’s experienced pushing) 0.10 0.05 [0.01, 0.21] NA NA NA NA NA 1.000 12092 10659
Additional Parameters
sigma 0.88 0.02 [0.85, 0.91] NA NA NA NA NA 1.000 59748 29706 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.84). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb +
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 4.00 [3.79,  4.23]         2.00      0.25
##    pressure_partner_cw 1.42 [1.37,  1.49]         1.19      0.70
##        pushing_self_cw 1.34 [1.29,  1.39]         1.16      0.75
##     pushing_partner_cw 1.16 [1.12,  1.20]         1.08      0.86
##     persuasion_self_cb 1.11 [1.08,  1.15]         1.05      0.90
##  persuasion_partner_cb 1.05 [1.03,  1.10]         1.03      0.95
##       pressure_self_cb 1.35 [1.30,  1.41]         1.16      0.74
##    pressure_partner_cb 1.17 [1.13,  1.22]         1.08      0.86
##        pushing_self_cb 4.57 [4.32,  4.83]         2.14      0.22
##       barriers_self_cb 4.88 [4.62,  5.16]         2.21      0.20
##   facilitators_self_cw 2.08 [1.99,  2.18]         1.44      0.48
##   facilitators_self_cb 2.28 [2.18,  2.40]         1.51      0.44
##                    day 1.34 [1.29,  1.39]         1.16      0.75
##  Tolerance 95% CI
##      [0.24, 0.26]
##      [0.67, 0.73]
##      [0.72, 0.78]
##      [0.83, 0.89]
##      [0.87, 0.93]
##      [0.91, 0.97]
##      [0.71, 0.77]
##      [0.82, 0.88]
##      [0.21, 0.23]
##      [0.19, 0.22]
##      [0.46, 0.50]
##      [0.42, 0.46]
##      [0.72, 0.78]
## 
## Moderate Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 8.03 [7.59,  8.51]         2.83      0.12
##  persuasion_partner_cw 9.56 [9.02, 10.13]         3.09      0.10
##     pushing_partner_cb 5.74 [5.43,  6.08]         2.40      0.17
##       barriers_self_cw 5.01 [4.74,  5.29]         2.24      0.20
##  Tolerance 95% CI
##      [0.12, 0.13]
##      [0.10, 0.11]
##      [0.16, 0.18]
##      [0.19, 0.21]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value = 0.52
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.00000000 0.00313786
## sample estimates:
## outlier frequency (expected: 0.000637860082304527 ) 
##                                         0.002057613
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 4.02*** 1.48 [ 2.00, 8.48] 1.000 [0.84, 1.20] 0.000 55.466 Very Strong Evidence 1.000 26171 29780
Intercept[2] 8.68*** 3.32 [ 4.21, 18.94] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 25869 28076
Intercept[3] 24.03*** 9.89 [ 11.12, 56.11] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 26031 27575
Intercept[4] 101.08*** 48.27 [ 41.38, 272.00] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 27952 29608
Intercept[5] 1719.09*** 1499.08 [386.52, 12530.15] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 43543 30931
Within-Person Effects
Daily individual’s experienced persuasion 0.74* 0.09 [ 0.57, 0.94] 0.994 [0.84, 1.20] 0.149 1.564 Weak Evidence 1.000 27412 27241
Daily partner’s experienced persuasion 1.02 0.13 [ 0.78, 1.30] 0.563 [0.84, 1.20] 0.828 0.060 Strong Evidence for Null 1.000 35721 28312
Daily individual’s experienced pressure 1.78* 0.40 [ 1.08, 2.80] 0.985 [0.84, 1.20] 0.049 1.077 Weak Evidence 1.000 27232 24546
Daily partner’s experienced pressure 1.24 0.40 [ 0.54, 2.41] 0.747 [0.84, 1.20] 0.322 0.103 Moderate Evidence for Null 1.000 22469 17310
Daily individual’s experienced pushing 1.32* 0.17 [ 1.03, 1.71] 0.985 [0.84, 1.20] 0.206 0.645 Weak Evidence for Null 1.000 40213 28324
Daily partner’s experienced pushing 0.90 0.14 [ 0.66, 1.22] 0.752 [0.84, 1.20] 0.666 0.072 Strong Evidence for Null 1.000 35218 26735
Day 1.42 0.68 [ 0.56, 3.57] 0.767 [0.84, 1.20] 0.227 0.068 Strong Evidence for Null 1.000 42741 30598
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.82 0.16 [ 0.56, 1.19] 0.857 [0.84, 1.20] 0.434 0.090 Strong Evidence for Null 1.000 54312 29838
Daily facilitators 0.83 0.15 [ 0.58, 1.18] 0.850 [0.84, 1.20] 0.468 0.091 Strong Evidence for Null 1.000 44095 30209
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.30 0.94 [ 0.31, 5.62] 0.639 [0.84, 1.20] 0.184 0.080 Strong Evidence for Null 1.000 20631 26275
Mean partner’s experienced persuasion 2.37 1.91 [ 0.52, 13.01] 0.866 [0.84, 1.20] 0.103 0.136 Moderate Evidence for Null 1.000 18794 23357
Mean individual’s experienced pressure 2.84 2.50 [ 0.50, 16.92] 0.883 [0.84, 1.20] 0.080 0.133 Moderate Evidence for Null 1.000 24507 28940
Mean partner’s experienced pressure 0.81 0.68 [ 0.14, 3.97] 0.598 [0.84, 1.20] 0.164 0.064 Strong Evidence for Null 1.000 19447 25763
Mean individual’s experienced pushing 0.76 0.78 [ 0.10, 6.28] 0.606 [0.84, 1.20] 0.133 0.074 Strong Evidence for Null 1.000 18595 22131
Mean partner’s experienced pushing 0.04** 0.05 [ 0.00, 0.45] 0.995 [0.84, 1.20] 0.004 2.235 Weak Evidence 1.000 20492 23212
Mean barriers 1.07 0.62 [ 0.34, 3.30] 0.544 [0.84, 1.20] 0.240 0.067 Strong Evidence for Null 1.000 28391 30272
Mean facilitators 0.92 0.32 [ 0.47, 1.86] 0.597 [0.84, 1.20] 0.390 0.081 Strong Evidence for Null 1.000 25625 26775
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.11 0.29 [0.63, 1.80] NA NA NA NA NA 1.000 11292 20153
sd(Daily individual’s experienced persuasion) 0.34 0.19 [0.03, 0.72] NA NA NA NA NA 1.000 5620 9730
sd(Daily partner’s experienced persuasion) 0.19 0.15 [0.01, 0.54] NA NA NA NA NA 1.000 13956 15401
sd(Daily individual’s experienced pressure) 0.47 0.32 [0.03, 1.27] NA NA NA NA NA 1.000 9034 12841
sd(Daily partner’s experienced pressure) 0.64 0.52 [0.03, 2.03] NA NA NA NA NA 1.000 8655 13191
sd(Daily individual’s experienced pushing) 0.22 0.17 [0.01, 0.63] NA NA NA NA NA 1.000 8627 13288
sd(Daily partner’s experienced pushing) 0.18 0.16 [0.01, 0.66] NA NA NA NA NA 1.000 16181 20013
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2]
##   (r = 0.79), b_Intercept[4] and b_Intercept[3] (r = 0.85),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.74). This
##   might lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="05_SensitivityBarriers-ONLY_DAYS_WITH_PLAN-_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    facilitators_self_cw + facilitators_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.10 [1.07, 1.14]         1.05      0.91
##  persuasion_partner_cw 1.22 [1.18, 1.27]         1.10      0.82
##       pressure_self_cw 1.03 [1.01, 1.09]         1.02      0.97
##    pressure_partner_cw 1.04 [1.02, 1.10]         1.02      0.96
##        pushing_self_cw 1.11 [1.08, 1.15]         1.05      0.90
##     pushing_partner_cw 1.22 [1.18, 1.27]         1.10      0.82
##     persuasion_self_cb 3.16 [3.00, 3.33]         1.78      0.32
##  persuasion_partner_cb 3.63 [3.44, 3.84]         1.91      0.28
##       pressure_self_cb 2.21 [2.11, 2.32]         1.49      0.45
##    pressure_partner_cb 2.09 [1.99, 2.19]         1.44      0.48
##        pushing_self_cb 2.50 [2.38, 2.63]         1.58      0.40
##     pushing_partner_cb 2.38 [2.27, 2.51]         1.54      0.42
##       barriers_self_cw 1.32 [1.27, 1.38]         1.15      0.76
##       barriers_self_cb 1.27 [1.23, 1.33]         1.13      0.78
##   facilitators_self_cw 1.26 [1.22, 1.32]         1.12      0.79
##   facilitators_self_cb 1.09 [1.06, 1.14]         1.05      0.92
##                    day 1.06 [1.03, 1.11]         1.03      0.95
##  Tolerance 95% CI
##      [0.87, 0.94]
##      [0.79, 0.85]
##      [0.92, 0.99]
##      [0.91, 0.98]
##      [0.87, 0.93]
##      [0.79, 0.85]
##      [0.30, 0.33]
##      [0.26, 0.29]
##      [0.43, 0.47]
##      [0.46, 0.50]
##      [0.38, 0.42]
##      [0.40, 0.44]
##      [0.73, 0.78]
##      [0.75, 0.81]
##      [0.76, 0.82]
##      [0.88, 0.94]
##      [0.90, 0.97]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         34%
##          beta         16%
##        cauchy         16%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 30 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value =
## 0.6217
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000520929 0.0114105115
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002057613
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.54 0.22 [0.24, 1.20] 0.933 [0.83, 1.20] 0.118 0.217 Moderate Evidence for Null 1.000 31068 31972
Within-Person Effects
Daily individual’s experienced persuasion 0.70* 0.11 [0.49, 0.94] 0.991 [0.83, 1.20] 0.118 1.426 Weak Evidence 1.000 28303 27783
Daily partner’s experienced persuasion 1.14 0.21 [0.80, 1.71] 0.771 [0.83, 1.20] 0.564 0.108 Moderate Evidence for Null 1.000 29693 27955
Daily individual’s experienced pressure 2.35* 1.10 [1.05, 9.09] 0.981 [0.83, 1.20] 0.040 1.158 Weak Evidence 1.000 17531 16805
Daily partner’s experienced pressure 1.83 1.31 [0.44, 15.11] 0.821 [0.83, 1.20] 0.141 0.240 Moderate Evidence for Null 1.000 23241 20120
Daily individual’s experienced pushing 1.44* 0.23 [1.07, 2.07] 0.992 [0.83, 1.20] 0.111 1.366 Weak Evidence 1.000 30892 30438
Daily partner’s experienced pushing 0.90 0.22 [0.54, 1.55] 0.671 [0.83, 1.20] 0.483 0.109 Moderate Evidence for Null 1.000 32781 26660
Day 1.31 0.71 [0.45, 3.79] 0.692 [0.83, 1.20] 0.230 0.067 Strong Evidence for Null 1.000 44502 29325
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.73 0.17 [0.47, 1.13] 0.919 [0.83, 1.20] 0.261 0.156 Moderate Evidence for Null 1.000 54569 30497
Daily facilitators 0.80 0.17 [0.53, 1.20] 0.856 [0.83, 1.20] 0.396 0.110 Moderate Evidence for Null 1.000 50848 31227
Daily wear time NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 3.52 3.70 [0.46, 31.79] 0.890 [0.83, 1.20] 0.067 0.224 Moderate Evidence for Null 1.000 21805 27694
Mean partner’s experienced persuasion 6.74 7.84 [0.78, 82.14] 0.958 [0.83, 1.20] 0.031 0.456 Weak Evidence for Null 1.000 20869 26194
Mean individual’s experienced pressure 56.85* 110.41 [1.59, 3857.92] 0.987 [0.83, 1.20] 0.007 1.534 Weak Evidence 1.000 19442 25397
Mean partner’s experienced pressure 0.47 1.03 [0.01, 28.93] 0.636 [0.83, 1.20] 0.063 0.171 Moderate Evidence for Null 1.000 15301 23026
Mean individual’s experienced pushing 0.51 0.85 [0.02, 15.63] 0.656 [0.83, 1.20] 0.077 0.125 Moderate Evidence for Null 1.000 16410 24931
Mean partner’s experienced pushing 0.02* 0.05 [0.00, 0.89] 0.979 [0.83, 1.20] 0.010 0.893 Weak Evidence for Null 1.000 17758 24064
Mean barriers 0.62 0.49 [0.12, 2.80] 0.733 [0.83, 1.20] 0.154 0.105 Moderate Evidence for Null 1.000 26724 28241
Mean facilitators 0.95 0.45 [0.38, 2.52] 0.541 [0.83, 1.20] 0.297 0.112 Moderate Evidence for Null 1.000 25700 28189
Mean wear time NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 2.02 0.43 [1.30, 3.05] NA NA NA NA NA 1.000 15283 23255
sd(Daily individual’s experienced persuasion) 0.47 0.22 [0.06, 0.94] NA NA NA NA NA 1.000 8335 9912
sd(Daily partner’s experienced persuasion) 0.42 0.24 [0.04, 1.01] NA NA NA NA NA 1.000 13773 13015
sd(Daily individual’s experienced pressure) 1.16 0.68 [0.10, 2.95] NA NA NA NA NA 1.001 8394 9744
sd(Daily partner’s experienced pressure) 1.69 1.06 [0.14, 4.25] NA NA NA NA NA 1.000 13251 12316
sd(Daily individual’s experienced pushing) 0.28 0.22 [0.02, 0.83] NA NA NA NA NA 1.001 11229 17180
sd(Daily partner’s experienced pushing) 0.45 0.33 [0.03, 1.29] NA NA NA NA NA 1.000 14798 17513
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length
## is not a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple
## of shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="05_SensitivityBarriers-ONLY_DAYS_WITH_PLAN-_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0     1.57      3.04     -0.4      5.4       5.77
##   Post.Prob Star
## 1      0.85     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate,
## stats_to_report = stats_to_report, : Coefficients were exponentiated.
## Double check if this was intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 2.46*** [ 1.48, 4.02] 49.83*** [42.49, 58.35] 122.48*** [108.02, 139.14] 3.80*** [ 3.57, 4.01] 0.54 [0.24, 1.20]
Within-Person Effects
Daily individual’s experienced persuasion 1.46*** [ 1.22, 1.78] 1.02 [ 0.97, 1.09] 1.02 [ 0.98, 1.06] 0.01 [-0.04, 0.06] 0.70* [0.49, 0.94]
Daily partner’s experienced persuasion 1.32** [ 1.10, 1.65] 1.02 [ 0.98, 1.07] 1.03 [ 0.99, 1.07] 0.01 [-0.04, 0.06] 1.14 [0.80, 1.71]
Daily individual’s experienced pressure 0.80 [ 0.47, 1.46] 0.88 [ 0.77, 1.01] 0.95 [ 0.86, 1.05] -0.04 [-0.18, 0.08] 2.35* [1.05, 9.09]
Daily partner’s experienced pressure 1.50 [ 0.81, 4.82] 0.95 [ 0.86, 1.04] 0.96 [ 0.90, 1.04] -0.06 [-0.19, 0.06] 1.83 [0.44, 15.11]
Daily individual’s experienced pushing 1.15 [ 0.86, 1.57] 0.97 [ 0.91, 1.03] 1.02 [ 0.97, 1.08] 0.04 [-0.04, 0.11] 1.44* [1.07, 2.07]
Daily partner’s experienced pushing 1.58** [ 1.17, 2.43] 0.97 [ 0.91, 1.02] 1.01 [ 0.96, 1.05] 0.11* [ 0.03, 0.18] 0.90 [0.54, 1.55]
Day 1.00 [ 0.65, 1.54] 1.00 [ 0.87, 1.14] 0.97 [ 0.88, 1.07] 0.24** [ 0.09, 0.39] 1.31 [0.45, 3.79]
Own actionplan NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.96 [ 0.80, 1.15] 1.21*** [ 1.13, 1.30] 1.01 [ 0.96, 1.05] -0.14*** [-0.20, -0.07] 0.73 [0.47, 1.13]
Daily facilitators 1.40*** [ 1.17, 1.68] 1.20*** [ 1.14, 1.26] 1.06** [ 1.02, 1.10] 0.12*** [ 0.06, 0.17] 0.80 [0.53, 1.20]
Daily wear time NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.86 [ 0.28, 2.49] 1.07 [ 0.78, 1.46] 1.02 [ 0.73, 1.42] 0.14 [-0.44, 0.73] 3.52 [0.46, 31.79]
Mean partner’s experienced persuasion 1.31 [ 0.43, 3.81] 0.92 [ 0.67, 1.26] 0.89 [ 0.64, 1.25] 0.45 [-0.13, 1.03] 6.74 [0.78, 82.14]
Mean individual’s experienced pressure 0.31 [ 0.07, 1.27] 0.64 [ 0.29, 1.40] 1.07 [ 0.74, 1.56] -0.16 [-0.76, 0.43] 56.85* [1.59, 3857.92]
Mean partner’s experienced pressure 0.25 [ 0.05, 1.05] 0.49 [ 0.23, 1.08] 1.02 [ 0.70, 1.45] -0.32 [-0.91, 0.27] 0.47 [0.01, 28.93]
Mean individual’s experienced pushing 1.13 [ 0.23, 5.71] 1.49 [ 0.83, 2.64] 1.11 [ 0.68, 1.80] 0.27 [-0.58, 1.11] 0.51 [0.02, 15.63]
Mean partner’s experienced pushing 1.42 [ 0.28, 7.17] 1.70 [ 0.94, 3.03] 1.29 [ 0.80, 2.08] -0.09 [-0.93, 0.76] 0.02* [0.00, 0.89]
Mean barriers 0.57* [ 0.37, 0.90] 0.99 [ 0.86, 1.15] 0.81*** [ 0.72, 0.89] -0.48*** [-0.64, -0.32] 0.62 [0.12, 2.80]
Mean facilitators 1.21 [ 0.89, 1.64] 1.11* [ 1.01, 1.22] 0.95 [ 0.89, 1.02] 0.29*** [ 0.18, 0.40] 0.95 [0.38, 2.52]
Mean wear time NA NA NA NA 1.00* [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 1.25 [0.94, 1.69] 0.30 [0.22, 0.40] 0.34 [0.26, 0.45] 0.62 [0.49, 0.81] 2.02 [1.30, 3.05]
sd(Daily individual’s experienced persuasion) 0.16 [0.01, 0.44] 0.11 [0.06, 0.17] 0.05 [0.02, 0.09] 0.04 [0.00, 0.12] 0.47 [0.06, 0.94]
sd(Daily partner’s experienced persuasion) 0.24 [0.02, 0.53] 0.06 [0.01, 0.12] 0.05 [0.00, 0.10] 0.04 [0.00, 0.11] 0.42 [0.04, 1.01]
sd(Daily individual’s experienced pressure) 0.43 [0.02, 1.79] 0.09 [0.00, 0.34] 0.07 [0.00, 0.23] 0.07 [0.00, 0.26] 1.16 [0.10, 2.95]
sd(Daily partner’s experienced pressure) 0.78 [0.04, 2.35] 0.05 [0.00, 0.19] 0.04 [0.00, 0.13] 0.09 [0.00, 0.29] 1.69 [0.14, 4.25]
sd(Daily individual’s experienced pushing) 0.45 [0.05, 0.96] 0.06 [0.00, 0.15] 0.09 [0.02, 0.16] 0.08 [0.00, 0.18] 0.28 [0.02, 0.83]
sd(Daily partner’s experienced pushing) 0.42 [0.03, 1.14] 0.05 [0.00, 0.13] 0.04 [0.00, 0.11] 0.10 [0.01, 0.21] 0.45 [0.03, 1.29]
Additional Parameters
sigma NA NA 0.66 [0.63, 0.69] 0.54 [0.52, 0.56] 0.88 [0.85, 0.91] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • rmarkdown (version 2.29; Allaire J et al., 2024)
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • rlang (version 1.1.4; Henry L, Wickham H, 2024)
  • tidybayes (version 3.0.7; Kay M, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • plotly (version 4.10.4; Sievert C, 2020)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()